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Abstract 

A general Monte Carlo simulation method of calculating the elastic constants of polydisperse 
hard-sphere colloidal crystal was developed. The elastic constants of a size polydisperse hard 
sphere fee crystal is calculated. The pressure and three elastic constants(Cn, C\i and C44) increase 
significantly with the polydispersity. It was also found from extrapolation that there is a mechanical 
terminal polydispersity above which a fee crystal will be mechanically unstable. 
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I. INTRODUCTION 



Elastic constants are among the most important physical quantities describing macro- 
scopic mechanical behaviors of a crystal. For hard sphere crystals, the elastic constants had 
been calculated by various groups with different methods including density functional theory 




0, j|, Monte Carlo simulations |g, H, S, S 
12] in the last two decades. The work of Jaric 
poisson ratio in the hard sphere crystal, which stimulated much interest in the investigation 
of the elasticity of this system. However, their predictions were proved incorrect by subse- 



10 and molecular dynamics simulations 
1] and Jones Q predicted a negative 



quent studies 



111 ] . The hard sphere system is the simplest in the systems which have 



pure repulsion interaction, thus the study of elasticity of hard sphere crystals is an excellent 
starting point to the studies of more complicated and realistic repulsive systems. Due to the 
simplicity of the hard sphere system, it often serves as a simple model system for testing new 
theoretical approaches. Hard sphere model is also an important model for a large class of 
colloid systems, understandings of colloid systems are greatly benefited from the extensive 
researches of the monodisperse hard sphere model. A more realistic model describing hard 
sphere colloidal systems is the size polydisperse hard sphere system. The size polydispersity 
of colloid particles is an intrinsic property of colloidal systems [3]. The polydispersity of 
colloidal hard spheres is characterized by the ratio of the standard deviation to the mean of 
the diameter. It has a remarkable influence on the thermodynamic and dynamic behaviors 



of the system 



H n, is, i7, Q : 
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221 ] . It is natural to expect that the elastic 



constants of a polydisperse hard sphere crystal differ from those of a monodisperse hard 
sphere crystal. However, it seems that the problem of elastic constants of polydisperse hard 
sphere crystals is not properly addressed in the literature to the best of our knowledge. 

In this work we propose a general Monte Carlo scheme to investigate the elasticity of 
a polydisperse hard sphere system, by which we calculated the elastic constants of a size 
polydisperse hard sphere crystal. We only consider face-center-cubic (fee) crystal structure 
in this study, as our previous calculations [23] showed that fee structure is still the most 
stable one for a size polydisperse hard sphere crystal as in the monodisperse case. When 
simulating a polydisperse crystal the semigrand ensemble is the best candidate 
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25]. 



In the ensemble the imposed physical quantity is the chemical potential difference Afi of 
particles of each kind to a reference kind. To obtain the chemical potential difference for a 
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prescribed size distribution p(er), one has to solve a functional inverse problem 



26 
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28|, 



which can be accomplished easily by the semigrand ensemble version of the nonequilibrium 
potential refinement (NEPR) method (SNEPR method) [23]. 

Simulation approaches for calculating elastic constants of a monodisperse hard sphere 
crystal are divided into two categories. One is the "fluctuation" method 7] where elastic 
constants are related to the thermal averages of the corresponding stress components. There 
are some difficulties to extend the method to the polydisperse system. In the present paper 
we employ another method, the so called "strain" method jf], 0, 10, 11, 12], where elastic 



constants can be obtained from the free energy-strain relation or its first derivative, the 
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29] 



stress-strain relation. In the simulation we used the extended ensemble method 
to determine the Helmholtz free energy of the crystal with different strain, which can also 
be obtained by thermodynamic integration method [30]. Then the elastic constants were 
extracted from the free energy-strain data. 

The paper is organized as follows. In Sec. II, we introduce the model and explain how 
the semigrand ensemble can be applied to calculate the elastic constants of polydisperse 
hard sphere crystals. Section III describes the simulation method employed in this work. 
The computational details and results are provided in Sec. IV. Finally, we present our 
conclusions in Sec. V. 



II. MODEL AND THEORY 

A. Semigrand canonical emsemble 

The semigrand canonical ensemble (SCE) is the most suitable ensemble for the simulation 
study of the elastic properties of a size polydisperse hard sphere crystal. In this ensemble the 
total number of particles and the distribution of particle sizes are fixed while the number of 
particles of each size is permitted to fluctuate. In the simulation study, the total number of 
particles in the system usually limited to a few hundreds to thousands, which are too small 
to distribute accurately according to a prescribed distribution of particle sizes. With the 
semigrand canonical ensemble, the distribution is realized through averages since the particle 
sizes are allowed to fluctuate. The grand canonical ensemble can also achieve the goal of dis- 
tribution realization but insertion and deletion of particles often require more computational 



3 



resources. For a size polydisperse system consisting of N particles with the composition dis- 
tribution p(a) in a constant volume V, the Helmholtz free energy F(N, V, {p(a)},T) of the 
system is given by 

F = -PV + N J p:{a)p{a)da, (1) 

where P is the pressure, a is the diameter of the particles, and p(cr) is the chemical potential 
of particles with diameter a. The semigrand canonical free energy(SCFE) Y(N, V, {p(a)}, T) 
is obtained from the Helmholtz free energy by the Legendre transformation, 

Y = F-N J (jjl((t) - p(a r ))p(a)da, (2) 

where a r is the diameter of an arbitrarily chosen reference component. The SCFE 
Y(N, V, {/i(cr)}, T) is a functional of p(cr) — p(cr r ). In semigrand canonical ensemble the 
thermodynamic variable to characterize the equilibrium system is /x(cr) — p(a r ) rather than 
the composition distribution p(a). 



The partition function 7 for this ensemble is 31] 



7 = M / " 1 Zs 



n A 3 (at 



N ~\ N 



X 



exp < (3 ^2(ji(<Ti) - n(<r r )) \ IJ da t , (3) 



i=l ) i=l 



here A(t7j) = /i/(27rmj/cT) 1//2 is the thermal wavelength of the ith. particle and is the 
canonical configuration integral for a given size configuration 

Z N = I ■■■ I e-wfldn. (4) 

Jri Jr N i=1 

By introducing the excess chemical potential relative to ideal gas p-exi^i) — ni&i) — 
kTln( NA y^ ), we can rewrite 7 in a more symmetric form 

7 = N i A 3N^ j / " / ZN X 6XP Y Y^ ex ^ ~ Vexi^r)) j *T<. (5) 

The partition function 7 is related to the thermodynamic potential Y through 

Y = -kT ln 7 (iV, V, T, p ex {a) - p ex {a r )). (6) 
In practical simulations the diameter of particles are discretized into a series of special 



particle sizes or divide the total size range of particles into many small bins 27|, |28( . As a 
result, the semigrand canonical partition function 7 becomes 



N 



^ max <-< max I I 

7 = mA 3N( a \ " ' Zn X exp \ P ^(Pesfa) - Vex(cr r )) I , (7) 

here <r m j n and a ma2: are, respectively, the minimum and maximum of particle sizes. 



B. Elastic constants of size polydisperse solid 



The elasticity of a solid is the property that the solid deforms in response to an external 
stress and return to its initial configuration when the stress is removed. Usually the de- 
formation is small, otherwise the deformation may become permanent. In the case of hard 
sphere solids, the situation is more complicated because an external stress is necessary to 
stabilize the solid, while the deformation of solid is induced by exerting an excess external 
stress. The deformation of a continuous solid can be described by the Lagrangian strain 
tensor 

^ 2 \dxj dxi dxi dxj ) ' 
here Ui is the ith component of the displacement field, and Xi is the ith component of the 
position of the displaced point in the solid, repeated indices are summed from 1 to 3. The 
Helmholtz free energy density / = F/V (where V is the volume of the undeformed solid) 
is a functional of the strain field. For a homogeneous deformation, namely the Lagrangian 
strain tensor 77^ is independent of the position, the free energy density becomes a function 
of the constant strain. The elastic constants can be defined in terms of the following Taylor 
expansion of the free energy density 

fivij) = /(o) + Tijio^ij + ^Cijmm + ■■■ , (9) 

where /(0) is the Helmholtz free energy density of the unstrained solid, and 7^(0) is the 
stress tensor of the unstrained solid which is necessary to stabilize the hard sphere crystal, 
in the case of an fee hard sphere solid, Tjj(0) = —p5ij where p is the hydrostatic pressure. 
The Cijki is the second-order elastic constants. 

It should be emphasized that in the case of a size polydisperse system, the composition 
distribution of the strained system must remain the same as that of the unstrained system, 
this is because 77^ and p{a) are two independent variables of the Helmholtz free energy. 
Thus the explicit expression of elastic constants is 

d 2 f(v) 



a 



ijkl 



(10) 

p(a);r)=0 



dVijdr] k i 

Here the subscripts p(cr) means that the distribution of particle sizes is fixed during the 
application of strain. In semigrand canonical ensemble, the Helmholtz free energy density 



5 



of the strained system with strain tensor r\ and composition distribution p{a) is 



N f 

f ^ N /"/ / a , ^^ 3NkT 
= y(V)+y / Wex(o, V) - Vex{(Tr, r]))p{a)da + - 



here y = Y/V. In the semigrand ensemble the composition distribution p(a) depends 
not only on the chemical potential difference but also on the strain rjij. Therefore, in 
order to keep the composition distribution unchanged, the chemical potential difference 
^ex{ a i v) ~ A*ea;( cr r) v) has to be adjusted for each strain 77^. That means at each given strain 
the chemical potential difference is recalculated. The evaluation of the chemical potential 
difference for a give n composition distribution and strain can easily be performed by the 



SNEPR method 



231 ]. explained in the following subsection. 



For the size polydisperse hard sphere fee crystal, there are only three independent second- 
order elastic constants Cmi, C1122 and Ci2i2- In the Voigt notation they can be written in a 
more compact format C\\ = Cuu, C12 = C1122 and C44 = 6*1212- It should be noted that the 



elastic constants defined in this way is not the one that is measurec 
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II 



directly in experiments, 



12|. 



though it is widely used in theoretical investigations |5|, [6|, [9|, 

For a cubic crystal (including fee structure) under isotropic pressure P the experimentally 
measured elastic constants, C T , are related to the above defined C by the following relations 

1 

Cn — C11 — P] C*i2 = C12 + P\ C44 = C44 — P- (12) 

It is only when P = that the two sets of elastic constants coincide. A detail discussion of 
the elastic constants under stress can be found in [321 ] . 



III. SIMULATION METHOD 



A. SNEPR method 



In order to fix the composition distribution of the polydisperse crystal, we extend the 
NEPR algorithm to the semigrand canonical ensemble. The algorithm can be used 
to find the chemical potential difference for a prescribed composition distribution, i.e. 
Ap ex (a) = A/i ex ({p(cr)}) at a given strain. Here, we only give a brief description of the 



method. A detailed presentation can be found in references 



23|, 



281 ] . For a given particle 




FIG. 1: (a), the solved excess chemical potential as function of the diameter of particles. Dotted 
line corresponds to the unstrained crystal, and solid line to the crystal with a contraction strain, 
(b), the solid line is the plot of the truncated Schultz function. The dotted line and dashed line 
are, respectively, the composition distribution for the unstrained and strained crystal, which are 
obtained from simulation by using the A(j, ex (a) plotted in (a). It is difficult to distinguish the three 
distributions. 

size distribution and strain, the chemical potential difference can be calculated by a Monte 
Carlo iteration procedure. Firstly, initial guess of the excess chemical potential is assigned 
to the implement of a semigrand canonical ensemble simulation, and then it is modified at 
every few MC steps according to the instant size distribution P ins (cr) as follows, 

A„» = A» ex (°) ~ 7, ( PmS pl~ {a ) ia) ) V - ( 13 ) 
Here 7; is a modification factor of the ith iteration. When the difference of the average size 
distribution P(a) and the given composition is less than a specified value £ 

P(a) - P(a) 



£ > max 



(14) 



P(a) 

one loop of the iteration is finished. The modification factor is then reduced, and the excess 
chemical potential of the last iteration is used as the initial input and start the next iteration. 
The iteration continues till the modification factor 7 reaches a very small value, typically 
1CT 5 , and the resulted excess chemical potential is then regarded as the solution of the 
problem. 

In order to check the validity of SNEPR method, we apply it to both unstrained and 
strained polydisperse hard sphere crystals. Figure 1(a) are the calculated excess chemical 



potential for the truncated Schultz distribution as function of the diameter of particles. The 
lower dotted line and upper solid line correspond to the unstrained crystal and the crystal 
with a contraction strain 0.004, respectively. In the contracted crystal it is more difficult 
to increase the particle size, so a higher excess chemical potential is necessary to fix the 
composition distribution. Figure 1(b) is the comparison of the given Schultz distribution 
and the distribution generated with the calculated excess chemical potential. In figure 1(b) 
we plotted three lines of the distribution, the solid line is the given distribution, dotted 
line is the distribution generated from the chemical potential difference for the unstrained 
crystal, and dashed line is the one for the strained crystal. The agreement among the three 
distributions is excellent and nearly undistinguishable from the figure. 



B. free energy difference calculation 

To determine the elastic constants one only needs to know the Helmholtz free energy 
difference between the unstrained system and the strained system with the same composi- 
tion distribution. The Helmholtz free energy difference can be decomposed into two parts 
according to equation ffTTl) 

A f(v) = A v(v) + ji^ex{^,ri) - n ex (a r ,T]))p(a)da^ , (15) 

where A denotes the excess quantities relative to the unstrained system. Correspondingly, 
the elastic constants becomes 

d 2 [Af( v )] 



C, 



ijkl 



drjijdrjki 



(16) 

p(o-);7?=0 



The second term on the right hand side of (|T5|) can be calculated when the fi ex (ai,r]) — 
Hex(o'r,v) f° r the fixed composition is determined. Our main task is thus to compute the 



is in- 



difference of semigrand free energy Ay(rj). To this end an extended ensemble [c 
troduced. The Lagrangian strain tensor r\ is regarded as an additional ensemble variable 
and different r] corresponds to different macroscopic state. The partition function of the 
extended ensemble is defined as 

Y}max 

F( N , T i P{°)) = ^2 ^ T ' /iei ^ cr ' ^ ~ aw(°v, v)), (17) 

r?=0 



s 



where r] max denotes the state with maximum strain and all macroscopic states 77 possess the 
same composition distribution p(a). From equation (jSJ) and (ITTj) the Ay{rj) becomes 

j(N, 0, T, p ex (a, 0) - fi ex (a r , 0)) 



Ay(r]) 



In 
In 
In 



T), T, /i ex (a, rf) - p ex (a r , 77)) 

7(Q)/r ' 
l/yfo)/r. 

Pr(0) 



LPr^j ' (18) 
where Pr(rj) is the probability that the system is in the macroscopic state rj. Therefore, 
knowledge of the macroscopic state probability distribution is sufficient to evaluate the elastic 
constants. The probability can be calculated from simulation by the flat histogram methods 
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361 ] . Here we demonstrate the implementation in the Wang-Landau scheme, other 



schemes can also be implemented. 

The extended ensemble Monte Carlo simulation involves three kind of moves, the first is 
the particle displacement, the second is the particle resizing and the third is the deformation 
of the simulation box which corresponding to the walk in the stain rj space. The first two 
moves are accepted or rejected in the usual Metropolis way, i.e. if no overlap between 
particles happens, the trial move (r, a) — > (r', a') is accepted with probability 



P aC c(r, a -> r', a 1 ) = min{l, exp(/3A// ea; (o-', 77) - f3Ap ex (a, 77))}. 



(19) 



Where Afi ex (a,r]) = fA ex (cr,r)) — fi ex (a r ,rj). The trial move in the r] space is treated with 
the Wang-Landau sampling in order to obtain the macroscopic state probability distribution 
Pr(rj). With a initial guess of the Pr(rj), the acceptance/rejection criteria for the simulation 
box deformation, 77 — ► 7/ is 



mm 




I Pr M (Yi) N p /9E w (A|*«(«T ) V)-AA* e .(<T ) i7)) 
Pr{rf) \ V ) e 



if no overlap of spheres 

otherwise 

(20) 

where rj only takes some discrete values 771, r)2---7) n - The chemical potential difference for each 
strain 77 were calculated with the SNEPR and stored before the simulation with extended 
ensemble. During the extended ensemble simulation the (unnormalized) Pr(rj) is updated 
by multiplying a modification factor / > 1 when a state of 77 is visited, a histogram of the 
distribution in 77 of the visited states is recorded to monitor the convergence of Pr (77). The 
relative probability distribution of 77 is obtained at the end of the simulation, which is then 



used to determine Ay(rj) from equation (|T8|) . Because the composition distribution p(er) is 
prescribed in advance, substituting p(cr), Ap ex (cr, 77) and Ay(rj) into f|T5|) the Helmholtz free 
energy density difference A/ (77) can then be determined. Finally, the elastic constants can 
be obtained from a polynomial fit to the free energy-strain data. 



IV. SIMULATION DETAILS AND RESULTS 



Most of the simulations are performed with a system of 256 size polydisperse hard spheres 
in a parallelepiped box, periodic boundary conditions are used in all three directions. The 
initial configuration is an ideal fee crystal (zero strain). The density of undeformed crystal 
taken in this simulation is p = 0.576. These value is chosen because there are some known 
results of the terminal polydispersity [ijj] and the elastic constants in the monodisperse case 
12] in literature. As the previous computations 10, 11] the elastic constants of fee 



crystal can be determined completely from three independent deformations. They are the 
contraction, contraction-expansion and shear deformation. The composition distribution of 
particles used in the simulation is the truncated Schultz function 

z+l 



z + l 



a 



a 



&min ^ ^ @maxi (^1) 



where a m i n and a max are respectively the minimum and the maximum diameters of hard 
sphere particles. W is the average diameter and z = 5~ 2 — 1 controls the width of the 
distribution. The criteria of the truncation is that the probability density at both ends of 
the distribution are almost equal. They are several times less than the maximum probability 
density, as shown in figure [T^b). Here the probability density at u m i n and cr max is 1/3 of the 
peak value. In the present paper the effect of the cutoff is not studied, the emphasis is on the 
effect of polydispersity to the elastic properties. The size polydisperse degree is defined by 

// _ — \2 

5 = v _ . In the simulation we consider a uniform discrete set of diameters of particles. 
When the number of discrete diameters is large enough, the diameter of particles tends to 
a continuous variable which can resemble the real polydisperse system. In this study 101 
different sizes of diameters were used. 

As mentioned above, the composition depends not only on the chemical potential dif- 
ference but also on the strain tensor. Therefore, the chemical potential differences are 
calculated before the extended ensemble simulation is performed. The results indicate that 



10 



the chemical potential difference increases with the magnitude of the strain, as plotted in 
Figure [T]( a). In the case of a contraction deformation, the increase of the chemical potential 
difference is more significant than the case of shear deformation. This is reasonable because 
the contraction deformation consumes more configuration space so that the turn to larger 
sizes is more difficult and requires a larger chemical potential. 

TABLE I: Elastic constants of fee hard sphere crystals and simulation parameters. All undeformed 
systems have the same density p = 0.576. Here, 5 is the size plydispersity, A the number of 
particles and P the pressure of the hard sphere crystal. The first two rows are the elastic constants 
of monodisperse hard sphere crystals obtained from reference |12j |. 



5 


A 


P 


C n 


Cl2 


C44 





256 




115.5(10) 


32.0(4) 


72.4(3.2) 





13292 




117.4(4.4) 


31.54(15) 


71.96(11) 





256 


14.54(2.5) 


115.1(46) 


32.7(16) 


72.0(35) 


0.03 


256 


15.49(2) 


125.1(47) 


43.7(14) 


77.1(45) 


0.03 


2048 


15.528(6) 


129.7(34) 


47.3(12) 


76.5(21) 


0.039 


256 


16.24(2) 


135.5(33) 


55.9(10) 


78.5(45) 


0.05 


256 


17.34(1) 


150.0(28) 


73.6(8) 


83.5(29) 



We performed simulations for four different size polydispersity, simulation parameters and 
results are given in Table [H The maximum polydispersity taken in our simulation is 5%, 

nnn 

because at higher S the crystal may be unstable [14], [l_5|, [38J. In order to check the system 



size dependence we performed simulations for a larger system with the particle number 
A" = 2048. The difference between the two systems are clearly less than or on the same 

ready large 



10 



12|]. We 



order of the statistical errors, which indicates that a system of 256 particles is a 
enough to get reasonable results, similar to the monodisperse hard sphere system 
also calculated the elastic constants of a monodisperse hard sphere crystal using the present 
method, which are in full agreement with the results of Pronk and Frenkel 12|. From Table 
U we see that the pressure increases with the increasing of the size polydispersity and the 
pressure with S = 0.05 is about 20 percent higher than the monodisperse system, which is 



consistent with the previous reports 



13, 



181 ] that the size polydispersity can induce a higher 



11 



osmotic pressure in the hard sphere colloidal crystal. A new phenomenon which has not 
been reported previously is that as the size polydispersity increases all elastic const ant s(Cn, 
C\i and C44) of a fee hard sphere crystal increase significantly. At 5 = 0.05, Ci 2 is 1.3 
times larger than that of the monodisperse system. Figure [2] shows the ratios of C\\ to 
C\2 and C44 to C*i 2 as a function of the polydispersity 5. It is interesting to note that for 
the polydispersity used in the simulation these ratios nicely follow a linear relation with the 
polydispersity. Experimentally, Phan et al have measured the high-frequency shear modulus 
for the hard sphere colloidal crystals [stJ. Their results are comparable to the static shear 
moduli of the fee crystal obtained in this study. Furthermore, our results indicate the 
static shear modulus also increases with the polydispersity. We think that the effect can be 
detected experimentally from a precise measurement. 

The monodisperse hard sphere system was regarded as a representative model in the 
description of many different aspects of hard sphere colloids. We expected that the poly- 
dispersity may give only a small correction to the monodisperse case. The large increase 
of the elastic constants with polydispersity indicates that the monodisperse model has its 
limitations in describing the real hard sphere colloids, especially in the elastic properties. 
The physics behind this large increase is still not clear, one general explanation is that the 
elastic constants are the second derivatives of the free energy, which should be more sen- 
sitive to the polydispersity than the free energy itself. A comprehensive understanding of 
this enhancement of elastic constants from polydispersity requires more extensive research 
which is beyond the scope of this paper. 

As is well known, the necessary requirements for the stability of a cubic crystal (including 
fee structure) are 

C£>0; C 4 r 4 >0; C£+2C£>0; (C^) 2 - {C^f > 0. (22) 

For the polydisperse hard sphere crystal under consideration, the first three conditions are 
obviously satisfied, the last one may be violated with increasing the polydispersity. Figure 
[3] depicts the change of — Cj 2 with the polydispersity. The value of Cf[ 1 — Cf 2 de- 
creases with increasing polydispersity, and tends to zero at 5 = 0.0807. To get this terminal 
polydispersity, we note that the simulation data follows a relation 

r T - r T - A - p aS+B 

°11 °12 — A e > 
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FIG. 2: The ratio of two elastic constants as a function of the size polydispersity 5. The triangles 
and squares indicate C\\/C\2 and C44/C12, respectively. The filled symbols denote the data ob- 
tained from 256-particle system, open symbols from 2048-particle system. The semi-filled symbols 
are from a 13292-particle system given by reference [jjj]. The lines are linear fits to the data. 

that means that In (A — Cf t + Cf 2 ) depends linearly on the polydispersity 5, as shown in 
figure [3J Using the fitted value of In A we can obtain the terminal polydispersity from 
extrapolation. This defines a mechanical terminal polydispersity (MTP) where the fee 
crystal becomes unstable. The strain related to the coefficient — Cj 2 is the following 
contraction-expansion deformation: 

x' = (l + e)x, y' = — 3,, z' = z. (23) 

That is to say, for 5 > 0.0807 the fee crystal is no longer stable under the deformation. 
The instability can also be described from the point of view of soft-mode. By solving the 
dispersion equation, one easily finds that there exists a transverse wave, propagates along 
the diagonal of the (001) crystal plane and polarized in x?/-plane, has the dispersion relation 
p m u) 2 = \(Cfi — Cf 2 )k 2 |39|], here p m is the mass density, u is the circular frequency and k 
the wave vector. Therefore, its frequency decreases substantially as the MTP is approached 
and this branch of the wave corresponds to a soft acoustical mode. 

lard sphere crystal, it is known that there is another terminal size-polydispersity 
38] above which the crystal will not be the most stable structure, and dis- 



For the 
Q, Q 
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order solid phase 



38| or solid-solid coexistence phase may become the most stable 
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5 

FIG. 3: The change of C£ - Cf 2 with the size polydispersity 5. The solid line is a linear fit to the 
data. The horizontal dotted line represents the function y = In (A) and A is a fitted constant(see 
text). At 5 = 0.0807 the solid line intersects with the dotted line. 

equilibrium state. We refer to the terminal polydispersity as thermodynamical terminal 
polydispersity (TTP). The MTP has to be not lower than the TTP, since above TTP the 
fee crystal exists in a metastable state, and has good mechanical behaviors. Therefore, the 
MTP obtained in this study gives an upper limit of the TTP. 

One important point of the semigrand canonical ensemble simulation of polydisperse 
hard sphere crystal is the swapping of particles. In this ensemble the particle sizes are 
allowed to fluctuate to get the required size distribution, with the fluctuation the particles 
of different sizes effectively exchange their spatial positions constantly, and with sufficient 
long simulation, the equilibrium state is realized during the simulation. The calculated 
elastic constants in this simulation is the equilibrium elastic constants, we may refer to 
them as ideal elastic constants. On the other hand, in real colloid crystals the particles can 
not exchange their positions simply because the free energy barrier is too high to be overcome 
in any reasonable time period. The particle arrangements are basically fixed by the process 
of crystal growth, which is not necessarily the equilibrium arrangement. It also noted that 
during the strain in the measurement of elastic constants the particles can only undergo 
small displacements, and it is not possible for particles to swap their positions. Based on 
this observation, a natural question is that wether the simulated elastic constants is the 
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same one in a measurement. The answer to this question is probably yes. The reason is that 
measurements are always performed with a crystal of macroscopic size, which contains much 
more colloid particles than the number of particles in the simulation study. If the sample is 
well relaxed, the self-averaging effect of macroscopic number of particles may compensate the 
effect of particle non-swap. We expect that the ideal elastic constant can be experimentally 
measured after sufficient equilibration. To test the possible deviation from the ideal elastic 
constants, we also performed simulations in the canonical ensemble. We randomly picked 
up several configurations from the particle sizes distribution, and for each configuration the 
particle size then fixed in the subsequent simulations of applying strains and extracting 
elastic constants. For a system of 2048 particles we find that the elastic constants with 
different configuration realizations of the same distribution give elastic constants within 
10% of difference. The average of the results from different realization of the configuration 
is listed in table [Til 

TABLE II: Comparison of the pressure and elastic constants obtained from the semigrand canonical 
ensemble simulation(upper row) and the canonical ensemble simulation (lower row). 





P 


Cii 


Cl2 


C44 


semigrand 


15.528 


129.7 


47.3 


76.5 


canonical 


15.591 


137.4 


49.9 


75.7 



V. CONCLUSION 

To conclude, the elastic constants of a fee polydisperse hard sphere crystal are simulated 
by the Monte Carlo method with a semigrand ensemble, the composition distribution is fixed 
in the simulation by the SNEPR method. The results show that both the pressure of the 
hard sphere solid and the three elastic constants increase with the size polydispersity 5. The 
method can be extended to the soft sphere system and the system with other polydisperse 
attributes in a straightforward manner. Our results also indicate that there is a MTP where 
the fee crystal is unstable, which provides us an upper limit of TTP. Above TTP we do not 
know which of the two structure (disorder solid and solid-solid coexistence phase) is the most 
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stable one. Determining the stable state from computer simulations requires more effort to 
accomplish. 
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